## Replication for Unilateral Powers, Public Opinion, and the Presidency by
## Andrew Reeves and Jon Rogowski Use this code to replicate
## Table 5 and Figure 3 in the article and Table 6 of the Online SI.


rm(list = ls())
## be sure to download this hack from the dataverse
source("plot.ci2.R")
library(plyr)
library(foreign)
library(car)
library(Zelig)
library(apsrtable)
load("rr_yougov.Rdata")
##
yougovData$pathwayfav <-
  recode(yougovData$immi_fav, "9=NA;8=NA;5=NA;1=4;2=3;3=2;4=1")
yougovData$ideology <-
  recode(yougovData$ideo5, "9=NA;8=NA;6=NA;1=5;2=4;3=3;4=2;5=1")
yougovData$approval <-
  recode(yougovData$obamaapp, "9=NA;8=NA;5=NA;1=4;2=3;3=2;4=1")
yougovData$eo_approval <-
  recode(yougovData$eo_approval, "9=NA;8=NA;5=NA;1=1;2=2;3=3;4=4")
yougovData$eo_guns <-
  recode(yougovData$eo_guns, "9=NA;8=NA;5=NA;1=1;2=1;3=0;4=0")
yougovData$eo_dream <-
  recode(yougovData$eo_dream, "9=NA;8=NA;5=NA;1=1;2=1;3=0;4=0")
yougovData$pid7 <-
  recode(yougovData$pid7, "99=NA;98=NA;8=NA;7=1;6=2;5=3;4=4;3=5;2=6;1=7")
## Table 5 in the main text
immigration.glm <-
  glm(
    eo_dream ~ approval + ideology + pathwayfav + eo_approval,
    data = yougovData,
    weights = weight,
    family = binomial(link = "logit")
  )
#summary(immigration.glm)
guns.glm <-
  glm(
    eo_guns ~ approval + ideology + eo_approval, data = yougovData, weights = weight, family = binomial(link = "logit")
  )
#summary(guns.glm)
eo.table <- apsrtable(
  guns.glm, immigration.glm,
  model.names = c("Research Gun Violence", "Defer Deportations"),
  Sweave = TRUE,
  coef.names = c(
    "Intercept", "Obama Approval", "Ideology", "Opposition to Executive Orders", "Approval of Pathway to Citizenship"
  )
)
## Table 6 from the online SI
immigration.apps.glm <-
  glm(
    eo_dream ~ approval + ideology + pathwayfav + eo_approval,
    data = subset(yougovData, yougovData$approval > 2),
    weights = weight,
    family = binomial(link = "logit")
  )
guns.apps.glm <- glm(
  eo_guns ~ approval + ideology + eo_approval,
  data = subset(yougovData, yougovData$approval > 2),
  weights = weight,
  family = binomial(link = "logit")
)
immigration.dapps.glm <-
  glm(
    eo_dream ~ approval + ideology + pathwayfav + eo_approval,
    data = subset(yougovData, yougovData$approval < 3),
    weights = weight,
    family = binomial(link = "logit")
  )
guns.dapps.glm <- glm(
  eo_guns ~ approval + ideology + eo_approval,
  data = subset(yougovData, yougovData$approval < 3),
  weights = weight,
  family = binomial(link = "logit")
)
eo.byapps.table <- apsrtable(
  guns.apps.glm, immigration.apps.glm,
  guns.dapps.glm, immigration.dapps.glm,
  model.names = c(
    "Research Gun Violence",
    "Defer Deportations",
    "Research Gun Violence",
    "Defer Deportations"
  ),
  Sweave = TRUE,
  stars = "default",
  coef.names = c(
    "Intercept", "Obama Approval", "Ideology", "Approval of Presidential Power", "Approval of Pathway to Citizenship"
  )
)
## Make Figure 3 in main text
yougovDataImm <-
  yougovData[,c("eo_dream", "approval", "ideology", "pathwayfav", "eo_approval", "weight")]
yougovDataImm <- na.omit(yougovDataImm)
z.out <-
  zelig(
    eo_dream ~ approval + ideology + pathwayfav + eo_approval,
    model = "logit",
    data = yougovDataImm,
    weights = yougovDataImm$weight
  )
x.out <-
  setx(z.out, approval = seq(from = 1, to = 4, by = .25))#, approval = 2.5)
s.out <- sim(z.out, x = x.out)

#pdf(file = "unilateral-paper/figures/yougov-plot.pdf", family = "Times")
layout(matrix(c(1,3,2,4), 2, 2, byrow = TRUE))
par(
  oma = c(1, 7.25, 2, 2),
  mar = c(3, 2, 0, 0),
  mgp = c(2, 1, 0),
  xpd = NA
)

plot.ci2(
  s.out,
  xlab = "Obama Approval",
  ci = c(95),
  ylim = c(0,1),
  legpos = c(200,200),
  ylab = "Probability of Support",
  main = "Immigration"
)
axis(
  side = 1, at = c(1,2,3,4), labels = c("","","","")
)
par(las = 1)
mtext("Obama\nApproval", 2, line = 4, cex = 1)
mtext("Strongly\n Approve", 1, at = 4, line = 1, cex = .5)
mtext(
  "Strongly\n Disapprove", 1, at = 1, line = 1, cex = .5
)
x.out <-
  setx(z.out, eo_approval = seq(from = 1, to = 4, by = .25))#, approval = 2.5)
s.out <- sim(z.out, x = x.out)

plot.ci2(
  s.out,
  xlab = "Support for Executive Orders",
  ci = c(95),
  ylim = c(0,1),
  legpos = c(200,200),
  ylab = ""
)
axis(
  side = 1, at = c(1,2,3,4), labels = c("","","","")
)
mtext("Presidential\nPower\nApproval", 2, line = 4, cex = 1)
mtext(
  "Strongly\n Disapprove", 1, at = 4, line = 1, cex = .5
)
mtext("Strongly\n Approve", 1, at = 1, line = 1, cex = .5)

### CDC Gun Control
yougovDataGuns <-
  yougovData[,c("eo_guns", "approval", "ideology", "eo_approval", "weight")]
yougovDataGuns <- na.omit(yougovDataGuns)
z.out <- zelig(
  eo_guns ~ approval + ideology + eo_approval,
  model = "logit",
  data = yougovDataGuns,
  weights = yougovDataGuns$weight
)
x.out <-
  setx(z.out, approval = seq(from = 1, to = 4, by = .25))#, approval = 2.5)
s.out <- sim(z.out, x = x.out)
plot.ci2(
  s.out,
  xlab = "Obama Approval",
  ci = c(95),
  ylim = c(0,1),
  legpos = c(200,200),
  ylab = "Probability of Support",
  main = "Gun Control"
)
axis(
  side = 1, at = c(1,2,3,4), labels = c("","","","")
)
mtext("Strongly\n Approve", 1, at = 4, line = 1, cex = .5)
mtext(
  "Strongly\n Disapprove", 1, at = 1, line = 1, cex = .5
)

x.out <-
  setx(z.out, eo_approval = seq(from = 1, to = 4, by = .25))#, approval = 2.5)
s.out <- sim(z.out, x = x.out)
plot.ci2(
  s.out,
  xlab = "Support for Presidential Power",
  ci = c(95),
  ylim = c(0,1),
  legpos = c(200,200),
  ylab = ""
)
axis(
  side = 1, at = c(1,2,3,4), labels = c("","","","")
)
mtext(
  "Strongly\n Disapprove", 1, at = 4, line = 1, cex = .5
)
mtext("Strongly\n Approve", 1, at = 1, line = 1, cex = .5)





####Predicted Probabilities: Immigration
yougovDataImm <-
  yougovData[,c("eo_dream", "approval", "ideology", "pathwayfav", "eo_approval", "weight")]
yougovDataImm <- na.omit(yougovDataImm)
z.out <-
  zelig(
    eo_dream ~ approval + ideology + pathwayfav + eo_approval,
    model = "logit",
    data = yougovDataImm,
    weights = yougovDataImm$weight
  )
set.seed(100)
x.out.low <- setx(z.out, approval = 1)#, approval = 2.5)
x.out.hi <- setx(z.out, approval = 3)#, approval = 2.5)
s.out <- sim(z.out, x = x.out.low, x1 = x.out.hi)
summary(s.out)

set.seed(100)
x.out.low <- setx(z.out, eo_approval = 1)#, approval = 2.5)
x.out.hi <- setx(z.out, eo_approval = 3)#, approval = 2.5)
s.out <- sim(z.out, x = x.out.low, x1 = x.out.hi)
summary(s.out)


####Predicted Probabilities: CDC
yougovDataGuns <-
  yougovData[,c("eo_guns", "approval", "ideology", "eo_approval", "weight")]
yougovDataGuns <- na.omit(yougovDataGuns)
z.out <- zelig(
  eo_guns ~ approval + ideology + eo_approval,
  model = "logit",
  data = yougovDataGuns,
  weights = yougovDataGuns$weight
)
set.seed(100)
x.out.low <- setx(z.out, approval = 1)#, approval = 2.5)
x.out.hi <- setx(z.out, approval = 3)#, approval = 2.5)
s.out <- sim(z.out, x = x.out.low, x1 = x.out.hi)
summary(s.out)

set.seed(100)
x.out.low <- setx(z.out, eo_approval = 1)#, approval = 2.5)
x.out.hi <- setx(z.out, eo_approval = 3)#, approval = 2.5)
s.out <- sim(z.out, x = x.out.low, x1 = x.out.hi)
summary(s.out)